Computing SNR for ECG Signals
Difficulty Level:
Tags pre-process☁ecg☁noise☁snr

A very important parameter to consider when analysing a signal is the Signal to Noise Ratio (SNR) - a metric that classifies objectively the quality of the acquisition, and like the name suggests, the relation between the intensity of the signal and the undesired noise in the acquired data, which is defined by: $ \\SNR= \frac{V_{pp}^{signal}}{V_{pp}^{noise}}\\ $, being $V_{pp}^{signal}$ and $V_{pp}^{noise}$ the peak-to-peak amplitude of the signal and noise segment, respectively.

To obtain this parameter, there is a big difference in the procedure when applying it to slow and rapid signals. Slow signals, as expected, have slow oscillations. To find the noise signal, you just have to subtract the filtered signal to the raw one. However, when it comes to rapid signals, such as the electrocardiography (ECG) signal, you would not have a correct noise signal this way.

In this Jupyter Notebook you will learn how to compute the SNR values for ECG signals, using samples that are available in our website.


1 - Importation of the required packages

In [1]:
# Importing specific functions from the numpy Python package
from numpy import ptp, zeros, mean

# biosignalsnotebooks python package
import biosignalsnotebooks as bsnb

# Package to calculate SNR with log
from math import log10

2 - Opening files and storing data in variables

If this step is challenging for you, try this Jupyter Notebooks on Opening a .txt file and Working with file headers .

2.2 - Opening ECG file

In [2]:
# ECG file
data_ecg, header_ecg = bsnb.load("../../signal_samples/supine_1min.txt", get_header=True)

# The loaded data is in the form of Python dictionary, where the first entry is the mac address 
# of the first device of the acquisition. This line, extracts that mac address.
mac_address = list(header_ecg.keys())[0]

# For this particular acquisition, the ECG data corresponds to the data in channel 1.
signal_ecg = data_ecg['CH1']

3 - ECG Analysis

Eletrocardiogram (ECG) is a signal that describes the electrical activity of the heart. It is one of the most known biological signals! This signal is repetitive and has many cycles per minute.

3.1 - Conversion of the raw signal to its physical meaning

The device used to acquire signals provides the raw data for post analyses. Thus, before starting the interpretation and processing, we will convert it to physical units, that in the case of voltage corresponds to Volts or milliVolts (V or mV) using the transfer function that is already implemented in the biosignalsnotebooks Python package. This transfer function may be accessed in the sensor datasheet .

In [3]:
# Get the sampling rate of the acquisition
sampling_rate_ecg = header_ecg['sampling rate']

# Get the resolution of the sensor
resolution = header_ecg['resolution'][0]

# Converting signal to mV (microVolts) unit
signal_mv = bsnb.raw_to_phy('ECG', 'biosignalsplux', signal_ecg, resolution, option='mV')

# Time axis generation
time_ecg= bsnb.generate_time(signal_mv, sampling_rate_ecg)

3.2 - Graphic Representation of raw signal

Notice the repetitive nature of the signal and also how fast it varies over time.

In [4]:
# Creating the graphic
bsnb.plot(time_ecg, signal_mv, title="Raw ECG", y_axis_label = "ECG Signal (mV)", x_axis_label = "Time (s)")

3.3 - Measuring Vpp signal : peak-to-peak amplitude of the signal component

Since the $V_{pp}^{signal}$ corresponds to the difference between the maximum and minimum values of the signal, we first must find them in order to proceed to the subsequent subtraction. Fortunately, Python and numpy offer functions to simplify this process.

In [5]:
# Finding the maximum and minimum values of the ECG signal
max_ecg = max(signal_mv)
min_ecg = min(signal_mv)

# Calculating the amplitude of the signal
vpp_signal_ecg = max_ecg - min_ecg

# Notice that this procedure is condensed in a single function in the numpy Python package:
vpp_signal_ecg = ptp(signal_mv)
In [6]:
print("Amplitude of ECG signal: {} mV".format(vpp_signal_ecg))
Amplitude of ECG signal: 1.8944408053621813 mV

The next plot is a visual representation of the procedure. The V pp value is given by the absolute difference between the lines that encompass the ECG signal.

In [7]:
#Another visual help to guide you
length_ecg = len(signal_mv)

max_line = zeros(length_ecg) + max_ecg
min_line = zeros(length_ecg) + min_ecg

bsnb.plot([time_ecg, time_ecg,time_ecg],[signal_mv, min_line, max_line], color=['#009EE3', '#E84D0E', '#E84D0E'],
          title="ECG Signal Amplitude", y_axis_label = "ECG Signal (mS)", x_axis_label = "Time (s)")

3.4 - Noise signal and measuring V pp noise : peak-to-peak amplitude of the noise component.

Here is an overview on how you should proceed to measure peak to peak amplitude for the noise component on repetitive and fast varying signals, specifically ECG signals:

  1. Slice the signal into several single-beat signals
  2. Compute the amplitude from all the areas where there is no presence of ECG signal and average all the obtained values to define the final Vpp noise value

3.4.1 - Finding R Peaks

The typical heartbeat in an ECG signal is composed by various waves that correspond to specific events in the heart. Specifically, there are the P, Q, R, S and T waves that occur in the heart during the contraction and relaxation of the atria, and contraction and relaxation of the ventricles. Being the R peak the most noticeable wave of each heartbeat, algorithms for heartbeat detection usually focus on finding these structures.
In this case, we will use the Pan-Tompkins algorithm, which allows to rapidly and accurately detect R peaks. For this, we will use the biosignalsnotebooks Python package.

In [8]:
time_r_peaks, amplitude_r_peaks = bsnb.detect_r_peaks(signal_mv, sampling_rate_ecg, time_units=True, plot_result= True)

Representation of noise segment in one beat

In order to explain the used method to calculate the SNR for ECG signals, we are giving a visual example using just one beat - we choose the fifth one. As previously mentioned, you should slice the noise segment in intervals to minimize errors when measuring noise signal amplitude. We broke it into 20 intervals. You do not have to worry about the code bellow as it is just a way to clarify the method to find the noise amplitude. Look at the graphic and notice:

  • One heartbeat of the ECG - meaningful values ;
  • The noise segment ;
  • Another heartbeat of the ECG - meaningful values .

In [9]:
# Time of the onset of the sixth R peak (remember that Python lists and arrays start at 0)
onset_sixth_hb = time_r_peaks[5]

# Time of the onset of the fifth R peak
onset_fifth_hb = time_r_peaks[4]

# The start of the noise corresponds to the interval between peaks. Through observation, the fifth heartbeat ends at
# around 0.5 s after its start, while the sixth starts at around 0.65 s after the start of the previous
time_start_noise = onset_sixth_hb + 0.5
time_end_noise = onset_sixth_hb + 0.65

# Then, we need to convert it to index to identify it in the signal. The values are cast to integers because all indexes are integers.
start_noise = int(time_start_noise * sampling_rate_ecg)
end_noise = int(time_end_noise * sampling_rate_ecg)

# Now we will identify the heartbeats. The procedure is analogous, and so we will do it in single lines.
# The parcels of 0.7 and 1.5 correspond to emprical values that need to be added to the R peaks in order
# to identify the onset of the heartbeats, once they do not start with the R peaks.
time_start = time_ecg[int((onset_fifth_hb + .7)*sampling_rate_ecg):start_noise]
beat_start = signal_mv[int((onset_fifth_hb+ .7)*sampling_rate_ecg):start_noise]

time_end = time_ecg[end_noise:int((onset_sixth_hb + 1.5)*sampling_rate_ecg)]
beat_end = signal_mv[end_noise:int((onset_sixth_hb + 1.5)*sampling_rate_ecg)]

# signal with noise values
time_noise = time_ecg[start_noise:end_noise]
beat_noise = signal_mv[start_noise:end_noise]
In [10]:
#Graphic 
bsnb.plot([time_start, time_end,time_noise],[beat_start, beat_end, beat_noise], color=['#009EE3', '#009EE3', '#11E868'], 
          title="ECG - One Beat", y_axis_label = "ECG Signal (mS)", x_axis_label = "Time (s)")

One of several intervals to measure amplitude of noise signal

In [11]:
#Another visual help to guide you
length_noise = len(beat_noise)

max_line = zeros(length_noise) + max(beat_noise)
min_line = zeros(length_noise) + min(beat_noise)

bsnb.plot([time_noise, time_noise,time_noise],[min_line, max_line, beat_noise], color=['#E84D0E', '#E84D0E', '#009EE3'],
          title="ECG - Noise Segment Amplitude", y_axis_label = "ECG Signal (mS)", x_axis_label = "Time (s)")

3.4.2 - Compute the average noise amplitude

To get a more precise value for the amplitude of noise, one should consider the noise present in the signal as a whole . However, once ECG signal is composed of high and low frequency components, it is not possible to just separate noise using either a low pass or a high pass filter. Thus, our approach is to identify every interval between heartbeats, that should correspond to flat lines with no signal. Thus, by getting the average value of the amplitude of those segments, we get an estimate of the amplitude of the noise in the signal.

In [12]:
vpp_noise_ecg = []

# For this task, we will follow the same procedure as shown before, but store the values in a list, so that we can then calculate the mean value.
for t in time_r_peaks:
    start = int((t + 0.5) * sampling_rate_ecg) # 0.5 - time between a peak and a flat 
    end = int((t + 0.65)* sampling_rate_ecg) # 0.65 time between a peak and the end of the flat
    interval = signal_mv[start:end]
    vpp = ptp(interval)
    vpp_noise_ecg.append(vpp)
    
vpp_noise_ecg = mean(vpp_noise_ecg)
In [13]:
print("Amplitude of noise signal: {} mV".format(vpp_noise_ecg))
Amplitude of noise signal: 0.04991307487582334 mV

3.5 - Computing SNR (dB)

Though the formula for the signal to noise ratio shown before is valid, it is more usual to represent this quantity in decibels (dB). This conversion is given by: $ \\SNR= 20 \times log_{10}{\frac{V_{pp}^{signal}}{V_{pp}^{noise}}}\\ $. Both quantities will be calculated.

In [14]:
snr_ecg = vpp_signal_ecg/vpp_noise_ecg

# The multiplication by 20 is because the signals are in the unit of (micro)Siemes
snr_ecg_db = 20 * log10(snr_ecg)

print("SNR for ECG signal: {}".format(snr_ecg))
print("SNR for ECG signal: {} dB".format(snr_ecg_db))
SNR for ECG signal: 37.95480062238765
SNR for ECG signal: 31.58533428813793 dB

SNR values are proportional to the quality of the acquisition. Thus, the higher the value, the better the acquisition, because the higher the influence of the signal relative to the noise. Specifically, the signal is higher than the noise if the SNR is higher than 1 (in dB higher than 0) and inversely, if it is below 1 (in dB below 0) the influence of the noise is higher than the influence of the signal and, thus, it might be impossible to recover the signal.

Though signal to noise ratio is important for every type of signal, the process for calculating it for the ECG signal is particular and was here demonstrated. We concluded that the noise influence was low on the overall signal, once the SNR is much higher than 1 (or 0 dB).

We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining Notebooks !

In [15]:
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()
.................... CSS Style Applied to Jupyter Notebook .........................
Out[15]: